Theoretical Population Biology

The height of ancestral trees, in other words the time to the most recent common ancestor in a genealogical tree, is an informative statistic in population genomics and used in contexts of genealogical dating, demographic history inference, and assessing signals of selection. Exploiting a general pathwise identity between the height of a tree and its branch lengths, we revisit Zeng et al.s mutation rate estimator to develop a model-free estimator for the product of the expected height of genealogical trees with the scaled mutation rate, both for single loci and genomic regions. This estimator is a linear function of the site frequency spectrum of a sample, i.e., of the distribution of allele counts at all segregating sites across the locus or genomic region and thus both observable and computationally cheap. We show that, under the infinite sites model of mutation, our estimator is unbiased and, for two genome-wide models of ancestries (sequential Markovian coalescents and common-pedigree ancestries), consistent. Furthermore, we show via simulation that our estimator performs smaller errors than averaging over reconstructed ancestral tree heights (extracted from reconstructed ancestral recombination graphs) when considering genome properties similar to human genomes across large genomic regions. We then revisit a publicly available genomic dataset of dogs and wolves and compare our method with extracting TMRCAs from ancestral recombination graph reconstructions.

The Yampolsky-Stoltzfus (YS) model provided the first formal demonstration that a bias in the introduction of new mutations can shape the direction of adaptive evolution, even when mutation rates are small relative to selection coefficients. Despite its foundational role, the model has been characterized almost entirely by simulation and closed-form solutions for its central quantity, the probability PA that a mutationally favored but less-fit genotype fixes ahead of a rarer, fitter alternative, have been available only in the limiting regimes of vanishing mutation supply and infinite population size. Here we derive closed-form approximations for PA across the intermediate regime of finite population size and appreciable mutation supply that lies between these two limits. Under strong selection between the competing mutants, we combine an establishment-probability argument with clonal interference theory to obtain an explicit expression for PA that accounts for both the ongoing input of the mutationally favored mutation and the frequency-dependent establishment probability of its competitor. Under weak selection between competing mutants, we derive a complementary expression based on the fixation of a mixed population under Kimura's classical formula. Both expressions agree closely with direct simulation of the underlying Wright-Fisher process, recover the known limiting cases, and remain consistent with one another over a substantial intermediate range of parameters. These results are expected to provide a basis for assessing how the strength of mutation bias affects adaptation in empirical data.

Single-cell sequencing provides a rich source of data that can be used to model and estimate the evolution of somatic cells within individuals. Many common models for modeling DNA evolution in phylogenetics are based on continuous-time Markov processes that operate on aligned sites. For haploid DNA, such as found in mitochondria, there are many well established models for these processes but since nuclear DNA is diploid, containing copies of both the maternal and paternal genome, it is natural to develop continuous-time Markov processes adapted to diploid DNA. One possible model is that both the maternal and paternal evolve independently according to haploid models so that, in principle, DNA comes with a phase. i.e. a label that tells if the site was derived from the maternal or paternal sequence. However, this information may not be available in practice and only an unphased diploid genotype is available. That is, an unphased model groups both of the possible maternal-paternal labels into a single state. This grouping process, called ``lumping'', is used by other methods as well and results in a coarser process. Care must be taken so that the coarser process can still be modeled as a continuous-time Markov chain; the lumped process of the original phased Markov chain may or may not be Markovian, depending on restrictions imposed on the generator matrix. Here, we explore the mathematical formulations of phased and unphased Markov processes, as well as a collection of other diploid DNA models from the literature.

Many pathogenic and microbial species undergo complex life cycles wherein they experience drastic bottlenecks during invasion into a host and rapid growth within hosts. The effects of such a complex population history on patterns of genetic variation is not currently understood. We model a pathogen population as a metapopulation where each deme represents a host and incorporate transmission dynamics as well as the effects of transmission bottlenecks. We employ a coalescent framework, obtaining analytical expressions for pairwise times to coalescence and genetic differentiation within and between hosts. We find that recurrent bottlenecks rescale the coalescent process within hosts to the size of the bottleneck, reducing within-host variation. In addition, the pairwise distribution of time to coalescence under our model deviates substantially from the Wright-Fisher process when the number of hosts (d) is of similar order to the size of the bottleneck (k), i.e., k [~] d and when m <<1/(2k) +1/(2d), where m is the migration rate. Otherwise, the pairwise process under our model converges in distribution to the Wright-Fisher process. Coalescent simulations show that under the same conditions the site frequency spectrum (SFS) under our model deviates from that under the Kingman coalescent and is drastically skewed towards singletons. Our results suggest that the inference of selection or demographic history from pathogen population genetics data may be systematically biased by the mechanics of pathogen reproduction, indicating a need for further theoretical work that can inform inference.

Population genetic dynamics are shaped by multiple evolutionary processes, including mutation, recombination, selection, and genetic drift. With the rapid growth of genomic data, modeling multilocus evolutionary dynamics and the resulting patterns of genetic variation has become increasingly important. Existing approaches often face challenges in jointly describing multiple evolutionary processes and modeling multilocus systems, motivating the development of more flexible and extensible frameworks. Here, we develop a variational framework for population genetic dynamics based on generalized gradient-flow theory, drawing on nonequilibrium thermodynamics. In our model, mutation, recombination, and selection are represented as modular variational components, with different life-cycle stages connected through a gamete-individual two-state system. Mutation and recombination act on the gamete state, selection acts on the individual state, and finite-population genetic drift is represented by a stochastic extension. Our model represents multilocus genetic variation directly in the full haplotype-frequency space, recovers classical mutation, recombination, and selection dynamics in the corresponding limits and provides a natural stochastic extension for finite populations. Numerical experiments and SLiM forward simulations show that our model captures the dynamics of allele frequencies, haplotype frequencies, and linkage disequilibrium in multilocus systems. The variational formulation opens avenues for future extensions to additional evolutionary processes and more complex multilocus systems, as well as for developing differentiable computational methods for gradient-based parameter inference and scalable genomic modeling.

Forward-in-time population genetic simulations are widely used in evolutionary analyses, but simulating large populations and long genomic regions remains computationally demanding. To reduce this cost, parameter rescaling is widely employed, in which the original evolutionary process is approximated by one with a smaller population size and fewer generations. Recently, several studies using the SLiM simulator have raised concerns about the accuracy of this rescaling approach. In this study, we show that many of the biases reported in these studies can be mitigated by using a different simulation algorithm. These results reveal that the accuracy of parameter rescaling depends on how well the simulation algorithm preserves diffusion-limit properties under rescaling.

Amongst measures of linkage disequilibrium, the correlation measure r2, due originally to Sewall Wright, has several advantages. Although basically a haploid statistic, it can readily be estimated using covariance and correlation from diploid data. Its population expectation can also be found from probability methods. We show that the probability of linked identity by descent, LIBD, directly estimates the expectation for r2 for Wright-Fisher models. The simplest such model for two loci is the haploid model, where each haplotype is chosen independently from haplotypes of the previous generation. This gives an equilibrium between recombination (recombination frequency c) and drift (population size N) where the expected value of r2 is approximately 1/(1+4Nc). For diploid models, the expectation is similar for low values of c, but contains a term in c2 and rises to double the haploid value for the limiting value of c = 0.5. An increase in r2 is also expected for the model of monogamous mating, as opposed to random mating for each offspring. The expectation for r2 is unaffected by the use of diploid rather than haploid data, an exception being the case of monogamous mating. Computer simulation shows close agreement of r2 with probability expectations, with minor discrepancies observed for unlinked loci. Using the probabilistic definition, r2 can be intuitively understood as a parameter shaping the effective sample size in genome-wide association studies.

Evolutionary rescue occurs when genetic adaptation prevents extinction of a declining population. Mutation and recombination generate adaptive variation, but when their rates at adaptive loci are too low, rescue may require evolution toward higher rates of genetic change. Here, we derive and compare analytical approximations for the probability of rescue through two alternative routes: mutation-rate and recombination-rate evolution. Our model comprises two resident genotypes from which the same adaptive haplotype can arise along either route, with rescue requiring the appearance of an unlinked modifier allele that increases mutation at, or recombination between, two fitness-determining loci. We show that resident genotype composition determines the relative success of the two routes. Mutation-rate evolution yields a higher rescue probability when one resident genotype is sufficiently more common because mutation can act on either genotype whereas recombination requires both. Otherwise, the two routes yield similar rescue probabilities. Even when either route makes rescue appreciably likely, however, the modifier alleles facilitating rescue are often lost before rescue is complete. This loss is driven by multiple factors: segregation can separate the modifier from the adaptive haplotype; the fitness cost of modifier expression raises extinction risk for modifier-bearing lineages; and continued modifier-induced mutation disrupts the adaptive haplotype after it forms. The latter effect can make recombination-rate modifiers more likely to persist even when mutation-rate evolution yields a higher rescue probability. Overall, these results show that both mutation- and recombination-rate evolution can induce rescue, but successful rescue need not entail persistence of the modifier allele that facilitated it.

We introduce a minimal individual based model for population intelligence that retains three structures absent from standard single agent descriptions: a manifold valued intelligence state, endogenous population turnover, and heritable lineage dependent reproductive fitness. The validated S^1 intelligence diffusion was extended with neutral logistic birth death dynamics. Each individual occupies an angle on the circle S^1; noisy mean field attraction produces collective synchronization; logistic birth death dynamics regulate population size; and reproductive rates can depend on lineage and position relative to designated safe or unsafe directions. This nested construction is intentionally small enough to admit an analytic equilibrium benchmark while remaining capable of selection, invasion, and safety transitions. Population size returned toward carrying capacity from (N(0)=200,800,1600), with tail means (775.7,800.3,805.1) for K=800. Across K=100 to 1600, the intelligence equilibrium error scaled as K^{-0.513} , consistent with the canonical K^{-1/2} fluctuation rate, while the extinction rate was zero in 120 runs over the simulated horizon. The mean order parameter remained close to the original theoretical (r*=0.7242), as neutral demography predicts. Directional reproductive selection shifts a safety score to +0.7447 or -0.7448 while synchronization remains approximately 0.75, demonstrating that coordination is not equivalent to safety. These are controlled numerical results, not asymptotic proofs.

The global potential of a chemical reaction network has many applications and is closely related to the stochastic detailed balance. However, many fundamental questions concerning stochastic detailed balance remain unresolved, such as whether it depends on the system volume and how to construct new systems that satisfy it. In this paper, we show that stochastic detailed balance may depend on the system volume. We therefore introduce four types of stochastic detailed balance according to their dependence on volume and rate constants, and systematically investigate the relationships among them. Our results distinguish detailed balance arising from particular choices of volume and parameters from that enforced by network structure, and identify conditions under which detailed balance at one volume extends to all volumes. We further obtain a class of networks satisfying stochastic detailed balance for every volume and every positive choice of rate constants, and construct new systems whose global potentials exhibit double-well structures.

Biological populations have been evolving in multi-dimensional trait spaces. Ecological interactions can induce the evolutionary diversification of a biological population into distinct multiple populations, which is called evolutionary branching in adaptive dynamics theory. If possible mutations are restricted to certain directions due to developmental, physiological, or physical constraints, the resulting evolutionary dynamics may be restricted to subspaces with fewer dimensions than the original trait spaces. In general, such constraints are also evolving slowly, due to rare and slight mutational deviations from the constraint subspaces. In this sense, constraints are not strict. However, the method for describing non-strict constraints that is globally consistent has yet to be developed. This study proposes a general and flexible method for constructing a non-strict constraint subspace from the corresponding strict constraint subspace. This method enables efficient analysis of adaptive evolution including evolutionary branching under non-strict constraints in the framework of adaptive dynamics theory. Applying this method to a resource competition model shows that the shape of the carrying capacity distribution can be treated as a non-strict constraint curve, and that the likelihood of evolutionary branching is affected not only by the shape of the constraint curve but also by its slow evolution.

Haldane's principle relates equilibrium mutation load to deleterious mutational input rather than to selection strength. We study this balance in a homeostatically regulated, damage-structured population where mutation occurs during reproduction, repair restores individuals toward lower-damage states, and nonlinear demographic feedback maintains constant population size. An exact demographic identity shows that load at every time equals excess reproductive effort weighted by mean relative fecundity. Without repair, the equilibrium branch with positive pristine-class frequency is exactly solvable, and along it the load is exactly the mutation probability times the realized per-capita reproductive turnover rate; when the terminal cap is removed, this is precisely the realized deleterious input, for every mutation probability below one. For rare mutation it reduces to the mutation probability times the pristine turnover rate, with the local fitness contrast canceling at leading order. With a finite terminal class the branch instead reaches a closure-dependent boundary collision. Repair modifies the rare-mutation balance by partitioning first-order removal from the first damaged class between selection and restoration, multiplying the leading load by the selective share of that removal. An exact binary solution extends these results to finite mutation pressure: the unrepaired interior equilibrium collides nonhyperbolically with the fully damaged boundary and converges algebraically at the threshold, whereas any positive repair rate restores a unique interior equilibrium for positive mutation; at the no-repair threshold, its displacement from the boundary scales as the square root of the repair rate as repair vanishes. These results separate the familiar competing-transition structure of mutation--selection balance from the specifically homeostatic effects of turnover-dependent mutation input and nonlinear demographic feedback.

Regular environmental cycles such as day and night or the lunar cycle are ubiquitous and drive many biological processes. Many organisms have evolved endogenous time-keeping mechanisms to anticipate these cycles to accurately time their behavior or synchronize physiological processes. The evolution of such synchronization to environmental cycles and their variability is not very well understood. As an example, we use the marine insect Clunio (Diptera, Chironomidae), which synchronize adult emergence and reproduction with the lunar cycle, and show different emergence patterns across populations. We use a model constructed with replicator-mutator dynamics to provide a novel outlook at the evolution of these emergence patterns. In particular, we investigate the effect of two traits of individual phenotypic variability of emergence, the spread and shape, on population-wide emergence patterns. We find that the semilunar and lunar emergence patterns are alternatively stable and occur when the phenotypic distribution has a narrow and is roughly flat. Arrhythmic emergence occurs when the phenotypic distribution is wide and nearly flat. The semi-lunar pattern can only be observed when the parameter of rhythmicity is changed. With our model, we introduce the concepts of frequency-dependent selection and individual phenotypic variability into the theory of chronobiology, which gives us insights into the evolution of rhythmicity in these organisms.

The Linkage Disequilibrium (LD) of unlinked loci can be used to estimate contemporary effective population size (Ne) of one to a few generations ago. In genomic datasets loci on different chromosomes are considered unlinked, but there are many more pairs of unlinked loci than there are independent pairs of chromosomes, resulting to confidence intervals (C.I.) being too narrow if the non-independence is not taken into account. Simulations were run to investigate the correlation structure among LD of unlinked loci, which can be expressed by the LD of loci along the same chromosomes, based on a discovery of a novel Random Probe LD estimator. We classify the correlation into two categories: overlapping of loci and disjoint pairs. The former is induced from the same locus being considered twice and is the stronger form of correlation. These correlations feed into {rho}, a parameter to quantify the degree of pseudo-replication in a dataset, and further a correction formula from which C.I. can be properly inferred. We demonstrate the use of our method via an analysis of genomic data from the malaria-transmitting Anopheles gambiae s.s mosquitoes. Apart from the point and C.I. estimates, we find that Var((r^2 ) ) is inflated by about 550 times due to pseudo-replication, highlighting the danger of not handling genetic correlation properly.

An important currency for individuals is their lifetime reproductive success (LRS), which is random simply due to demographic stochasticity. That randomness determines extinction probability. However, the distribution of LRS is also significantly affected by environmental variation. Previously we have shown (for a random environment that follows a Markov chain) how LRS is affected by an individual's birth environment. But our previous analysis severs the temporal linkage between a parent's birth environment and the environments into which its offspring are born. Here, we show how to compute the exact joint probability distribution of LRS for lineages spanning multiple environmental states (assuming a Markovian environment). From this joint distribution, we derive exact lineage extinction probabilities that fully incorporate demographic stochasticity, environmental frequency, and temporal autocorrelation. Applying our framework to Pacific Chinook salmon (semelparous with extreme early mortality) and European roe deer (iteroparous with delayed maturity), we demonstrate that initial birth states shape lineage fate. For salmon, the initial environment permanently separates trajectories; a poor birth state leads to near-certain extinction regardless of subsequent environmental shifts. For roe deer, we discover a counterintuitive dynamic where highly persistent poor conditions can rescue highly vulnerable individuals by expanding the right tail of reproduction. Furthermore, our exact calculations reveal that aggregated one-dimensional LRS distributions homogenize the reproductive landscape, overestimating extinction risks for lineages originating in poor environments and underestimating them for those in good environments. Accurately predicting evolutionary viability and the establishment of advantageous mutations requires preserving the environmental covariance. As global climate change amplifies environmental volatility, utilizing exact joint demographic models is critical for assessing true extinction risks.

Data augmentation is now a standard device across capture--recapture and occupancy analysis: adding a fixed number M of all-zero encounter histories replaces a model of unknown dimension with one of fixed dimension. Although M is often treated as a computational tuning choice, it also specifies a finite superpopulation and hence a binomial support constraint on the number of undetected individuals. In a Bayesian implementation that constraint appears as an induced prior; in a likelihood implementation it is the same finite-support assumption reached by another route. That the Bernoulli specification for the inclusion indicators induces a binomial prior on abundance is established (Schofield & Barker 2014); our concern is what that choice costs in estimated uncertainty. We develop the argument using a simple closed-population abundance estimation problem. We show that augmented occupancy and Huggins conditional-likelihood analyses give numerically identical point estimates of N once M is sufficiently large. Their uncertainty estimates, however, need not agree. We distinguish two sources of discrepancy. First, when M is small relative to the number of undetected individuals, the finite binomial ceiling truncates the likelihood or posterior and suppresses uncertainty. Second, once that ceiling no longer binds, Taylor-series (Delta-method) approximations still understate variance, because the quantity of interest is a strongly non-linear function of the estimated parameters and local linearization does not reproduce its curvature. Gauss-Hermite quadrature on the unconstrained logit scale recovers much of the shortfall and approaches the MCMC posterior benchmark, though a small residual remains that does not close as M grows, reflecting the distinction between asymptotic likelihood theory and finite-sample Bayesian inference. Neither mechanism is peculiar to abundance estimation: the first follows from the augmented representation itself, the second from any derived quantity that is a non-linear function of estimated parameters. We close with framework-specific guidance for choosing M.

In studies of assortative mating, similarity between variables measured in parents is often quantified using correlation. The order of the parents within any given pair can be arbitrary in these applications, but common correlation estimators are not robust to reordering within pairs. These unordered variable pairs are exchangeable, since the joint distributions of both orders are equal, and a given order is biased if the one variable has a lower expectation than the other. In this work, we characterize the effect of order bias on Pearson correlation estimates assuming exchangeable variables, and develop a new unbiased estimator, CorSym, that does not depend on order within each pair. Exchangeable variables have equal marginal distributions for both variables, a property accounted for by CorSym. In contrast, standard correlation estimators assume the two variables have different distributions, so biased orders skew the underlying mean, variance and covariance estimates. We show, through theory and simulations, how order bias often results in upwardly biased Pearson correlation estimates. Simulations confirm CorSym is unbiased, and validate its estimated confidence intervals. Using real admixed trios (parents and a child) from 1000 Genomes, we first demonstrate that the global ancestry of fathers and mothers are consistent with exchangeability, using both Kolmogorov-Smirnov tests and a Binomial test for order bias. However, ANCESTOR, which estimates parental global ancestry from a child's local ancestry, produces significant order biases in its output that result in substantial Pearson biases, which CorSym overcomes. Compared to ancestry proportions calculated directly on the parents, ANCESTOR also overestimates parent ancestry divergence and experiences another estimation artifact. Overall, CorSym solves an important estimation bias likely to be encountered in the study of assortative mating, providing unbiased and deterministic estimates that do not depend on the arbitrary order of the data.

Many questions in population genetics are concerned with reconstructing evolutionary history through time, such as inferring how population structure has changed throughout the past. Yet, many existing approaches have only an implicit temporal component, using quantities such as allele frequency or haplotype length as rough proxies for age. Recent advances in the inference of Ancestral Recombination Graphs (ARGs) have made it possible to estimate the entire sequence of local genealogies along the genome. These genealogies explicitly encode how samples are related to each other at different time points in the past, enabling the inference of how population structure has changed over time. To this end, recent work has used ARGs to define time-stratified versions of widely-used population genetics summary statistics in an attempt to capture the population structure present within a particular time window. Here, we show that naive approaches result in statistics that cannot be interpreted solely in terms of the population structure present within the time window they are targeting. To address this problem, we introduce a framework of coalescent-based time-stratified statistics, which use coalescence probabilities to partition classical summary statistics into interval-specific contributions. Using coalescent simulations, we demonstrate that these statistics accurately isolate population structure at different temporal depths and avoid spurious signals. Our results highlight the necessity of integrating coalescent theory into ARG-based temporal analyses and provide a principled and practical foundation for studying the dynamics of population structure through time.

We consider polygenic divergent selection in a mainland-island model, where our aim is to understand how patterns of genetic variation along the genome reflect the genetic architecture of postzygotic reproductive isolation. We derive a new expression for the effective migration rate () at both neutral and divergently selected loci (i.e. barrier loci), and develop a numerical approach to co-predict and allele frequencies at barrier loci. Using this , we predict neutral coalescence times along the genome, and show how our results for the mainland-island model can be used to obtain coalescence time predictions for nonequilibrium demographic models. We validate our approach extensively using forward-in-time individual-based simulations and show that our -based approximation remains accurate across a broad parameter range. We study both the case of a single barrier locus and a pair of barrier loci in detail to identify the conditions under which an -based approach performs well. We use our new predictions to examine how the genetic architecture of divergent selection shapes barriers to gene flow, quantifying hitchhiking effects along the genome and evaluating the effects of polygenicity on reproductive isolation. We discuss the implications for mapping barriers to gene flow using population genomic data, both for model-based inference and genome scan approaches based on summary statistics.

Functional protein sequence space is often described either by the sparsity of functional sequences or by the connectivity of neutral networks. These quantities characterize different properties: density measures the fraction of all sequences that are functional, whereas connectivity describes relationships among functional sequences. Neither alone determines how much of sequence space lies close to function. The coverage function measures the proportion of sequence space within a prescribed Hamming distance of functionality and therefore provides a direct geometric measure of local accessibility. Here we connect this geometric framework to characteristic discovery times using an explicitly heuristic model of stochastic exploration. The explored region is represented by an effective Hamming ball, and radial displacement is modelled as an outward-biased substitution process. A numerical illustration, calibrated to an empirical human germline mutation rate, shows how a coverage radius can be converted into a timescale. We then extend the framework to parallel search by defining an overlap-adjusted effective number of trajectories. The purpose is not to predict exact evolutionary waiting times, but to separate the geometric distribution of function from the dynamics by which sequence space is explored.